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ABSTRACT 

We investigate stationary axisymmctric configurations of magnetized stars in the framework of gen- 
eral relativistic ideal magnetohydrodynamics. Our relativistic stellar model incorporates a toroidal 
magnetic field and meridional flow in addition to a poloidal magnetic field for the first time. The 
magnetic field and meridional flow are treated as perturbations, but no other approximation is made. 
We find that the stellar shape can be prolate rather than oblate when a toroidal field exists. We 
also find that, for fixed baryonic mass and total magnetic helicity, more spherical the star is, lower 
the energy it has. Further, we find two new types of the frame dragging effect which differ from the 
^ standard one in a rotating star or Kerr geometry. They may violate the reflection symmetry about 

the equatorial plane. 

Subject headings: gravitation — MHD — relativity — stars: interiors — stars: magnetic fields - 
jy^ ' stars: neutron 
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1. INTRODUCTION 



Equilibrium configurations of stars are fundamental to study their dynamics, such as the precession, oscillation and 
CI gravitational wave emission (Ioka 2001; Bonazzola & Gourgoulhon 1996; Cutler 2002; Ioka & Taniguchi 2000). Among 
■ the main factors that affect the stellar structure, the magnetic stress may be very important especially for neutron 
J""^ ' stars (Bocquet et al. 1995; Bonazzola & Gourgoulhon 1996; Ioka 2001; Konno, Obata, & Kojima 1999). While most 
neutron stars have magnetic fields of ~ I0 12 -10 13 G, there is growing evidence for the existence of super-magnetized 



neutron stars with ~ 10 14 -10 15 G, the so-called magnetars, and their birthrate is estimated to be high, > 10% of 
all neutron stars (Kouveliotou et al. 1998; Duncan & Thompson 1992; Mereghetti et al. 2002; Thompson 2001). 
Larger magnetic fields > 10 16 G may be generated by the helical dynamo inside a new born neutron star (Duncan 
& Thompson 1992; Thompson & Duncan 1993). Even the maximum field strength allowed by the virial theorem 
Oh r*j 10 18 G could be achieved if the central engine of gamma-ray bursts are magnetars (Usov 1992; Nakamura 1998; 
q Kluzniak & Ruderman 1998; Wheeler, Hoflich, & Wang 2000). In a previous paper (Ioka & Sasaki 2003), we presented 
i_i ' a formalism to study the equilibrium configurations of magnetars. In this paper, based on this formalism, we calculate 
^ \ the actual configurations of magnetars. 

There has been some work on stationary axisymmetric configurations of magnetized stars in the framework of general 
relativistic magnetohydrodynamics (MHD) , but allowing the existence of only a poloidal magnetic field (Bonazzola et 
al. 1993; Bocquet et al. 1995; Konno, Obata, & Kojima 1999; Cardall, Prakash, & Lattimer 2001). The reason is that 
the existence of only a poloidal field is compatible with the circularity of the spacetime (Gourgoulhon & Bonazzola 
' 1993; Oron 2002), and once the circularity is assumed, the spacetime metric becomes substantially simple. In a circular 
| spacetime, there exists a family of two-surfaces everywhere orthogonal to the plane defined by the two Killing vectors 
associated with stationarity rf L = (d/dt)^ and axisymmetry £ M = (d/dip)^ (Papapetrou 1966; Carter 1969, 1973). 
Thus one may choose the coordinates (x^) = (t, x 1 , x 2 , ipj such that the metric components goi, go2, 331 and 532 are 
identically zero. As a consequence, the problem is simplified dramatically. However non-negligible toroidal magnetic 
fields are likely to exist in nature. In addition, a meridional flow may also exist in the interior of a neutron star, which 
also violates the circularity of the spacetime (Gourgoulhon & Bonazzola 1993; Oron 2002). Thus, we have to consider 
noncircular spacetimes. 

The problem to obtain an equilibrium configuration of a magnetized star can be separated into two parts. The first 
part is the matter and electromagnetic field equations in a given spacetime geometry (e.g., Zanotti & Rezzolla 2002; 
Rezzolla, Ahmedov, & Miller 2001a, b). The second part is the Einstein equations which determine the spacetime 
geometry under a given configuration of matter and electromagnetic fields. In Ioka & Sasaki (2003), we formulated 
the first problem, i.e., the equations of motion for the matter and electromagnetic fields in a curved spacetime. We 
reduced basic equations to a single differential equation, the so-called Grad-Shafranov (GS) equation for the magnetic 
flux function in a noncircular (i.e., the most general) stationary axisymmetric spacetime. 

In this paper, we develop a relativistic stellar model in which both a toroidal magnetic field and meridional flow 
are incorporated, in addition to a poloidal magnetic field, in a self-consistent way for the first time. We solve the GS 
equation and the perturbed Einstein equations simultaneously to determine the effects of toroidal fields and meridional 
flows on the spacetime geometry and stellar structure. We assume that magnetic fields are weak compared with gravity. 
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This assumption is valid as long as the magnetic field strength is smaller than the maximum value allowed by the 
virial theorem ~ 10 18 G (Bocquet et al. 1995; Bonazzola & Gourgoulhon 1996). We therefore treat the magnetic field 
as a small perturbation on an already-known non-magnetized, non-rotating configuration. Our approach is similar to 
that developed for slowly rotating stars (Chandrasekhar 1933; Hartle 1967; Hartle & Thorne 1968; Chandrasekhar & 
Miller 1974), in which the perturbation parameter is the angular velocity, whereas it is the magnetic flux function in 
our case. 

This paper is organized as follows. In §[21 we briefly review general relativistic ideal MHD in stationary axisymmetric 
spacetimes and recapitulate the GS equation obtained in Ioka & Sasaki (2003). In §|3 we take a weak magnetic field 
limit of the GS equation. This makes it possible to solve the GS equation by separation of variables. In §0J we derive 
the perturbed Einstein equations for the metric perturbations. In §0 we numerically solve the GS equation and the 
perturbed Einstein equations to obtain the magnetic field and fluid flow structure, the mass shift due to the magnetic 
field, the deformation of stars, and the frame dragging effects. Finally, we summarize our results in §EJ 

We use the units c = G = ks = 1- Greek indices (fj,, v, a, (3, ■ • •) run from to 3, small Latin indices (i, j, k, • ■ •) 
from 1 to 3, and capital Latin indices (A, B, C, ■ ■ •) from 1 to 2, where x° — t and x 3 = (p. The signature of the 
4-metric is (— , +, +, +). 

2. GENERAL RELATIVISTIC IDEAL MHD IN STATIONARY AXISYMMETRIC SPACETIMES 

The ideal MHD equations for a stationary axisymmetric system can be reduced to a single equation, the relativistic 
Grad-Shafranov (GS) equation (Ioka & Sasaki 2003; Punsly 2001). It is a second-order, nonlinear partial differential 
equation for a quantity called the flux function VP. The flux function "J is such that it is constant over each surface 
generated by rotating the magnetic field lines (or equivalently the flow lines) about the axis of symmetry and the GS 
equation determines the transfield equilibrium. Any physical quantities can be calculated from the solution ^> of the 
GS equation. The relativistic GS equation in a noncircular (i.e., the most general) stationary axisymmetric spacetime 
was derived in Ioka & Sasaki (2003). In this section, we briefly review the derivation. The weak magnetic field limit 
of the GS equation is given in § 

2.1. General relativistic MHD equations 

The basic equations for general relativistic ideal MHD are as follows (Lichnerowicz 1967; Novikov & Thorne 1973; 
Bekenstein & Oron 1978). Baryons are conserved, 

(p^)^ = 0, (1) 

where p is the rest mass density (i.e., the baryon mass times the baryon number density) and is the fluid 4- velocity 
with u^u^ = — 1. The electromagnetic field is governed by the Maxwell equations, 

F\jw,a] = 0j (2) 

F» v . v = A-kJ», (3) 

where and J M are the field strength tensor and the electric current 4- vector, respectively. Equation J2J) implies 
the existence of a vector potential A^; F M „ = A v>fl — A^, v . The electric and magnetic fields in the fluid rest frame are 
defined as 

E^F^u", (4) 
B^-^e^^F^, (5) 

where e M „ Q/ g is the Levi-Civita antisymmetric unit tensor with £0123 = \J— 9- In the ideal MHD, we assume the perfect 
conductivity, so that 

E lt =F ltv u v =0. (6) 

The equations of motion for the fluid are given by T MI/ ;t/ = 0, where is the total energy momentum tensor of the 
fluid and the electromagnetic field, 



= {p + pe+ y)u»u v + pg^ + ±- 

47T 



//"//" • 1//"' ) If 



(7) 



Here e and p are the internal energy per unit mass and pressure, respectively, B 2 := B^B V and we have used the 
perfect conductivity condition (jSJ). Assuming local thermodynamic equilibrium, the first law of thermodynamics is 
given by 

de = -pd(-) + TalS, (8) 



where S and T are the entropy per unit mass and the temperature. Finally we supply the equation of state, 

p=p(p,S). (9) 



Rclativistic stars with magnetic fields and meridional flow 



3 



and 



2.2. Conservation laws in the axisymmetric stationary case 

We take = (d/dt)^ and £ p = (d/dtp) 1 * so that x° = t and x 3 = <p are the time and azimuthal coordinates 
associated with the Killing vectors rf and £ M , respectively. Thus all physical quantities are independent of t and ip. 

In a stationary axisymmetric MHD system, there exist five conserved quantities along each flow line constructed 
from the energy-momentum tensor (Bekenstein & Oron 1978, 1979; Ioka & Sasaki 2003): D, L, 0, C and S. By 
exploiting the gauge freedom to make A^^rf — — and A^ v t^ = A^, 3 = 0, we can show that the magnetic 
potential W := A^ = A 3 as well as the electric potential $ := A^rf = A are constant along each flow line, i.e., 
tt i = = 0. Henceforth we label the flow line by 'J, which we will refer to as the flux function. The \& = const, 

surfaces are called the flux surfaces, which are generated by rotating the magnetic field lines (or the flow lines) about 
the axis of symmetry. Then one can show that 

^03-0, (10) 

F 0A = SIFa3, (11) 

F 3 i=-*,i = C^g~pu 2 , (12) 

F 23 = ^,2 = C^pu\ (13) 

F 12 = CV=gp(u 3 -nu°), (14) 

where and C(\I/) are conserved along each flow line and hence are functions of the flux function 'J. It may be 

useful to rewrite the above equations as 

B» = -Cp[{u + fiu3K + rf + fie]- (15) 

In addition we can show that E(^f), L(\&) and £>(*) are also conserved along each flow line where 

- D = f i(u + nu 3 ), (16) 

- E = (" + €p) m + C{U ° + QU3) ^ (17) 
L=^+^^«3 + C(uo + n«3)^, (18) 

/i = l + e+-, (19) 
P 

is the enthalpy per unit mass. These conserved quantities are not mutually independent but there is a relation among 
them. 

D = E-SIL. (20) 

Except for the entropy per unit mass S(^f), there are no perfectly relevant physical interpretations of these quantities. 
Nevertheless, by considering several limiting cases, we may associate them with terms that describe their qualitative 
nature. We may call 23 the fluid energy per unit mass, E(\&) the total energy per unit mass, L(^f) the total 
angular momentum per unit mass, tt(*f?) the angular velocity, and C(%?) the magnetic field strength relative to the 
magnitude of meridional flow. Since these conserved quantities are essentially the first integrals of the equations of 
motion, specification of these functions characterizes the configuration of the electromagnetic field and fluid flow. 

2.3. (2 + 1) + 1 formalism 

To describe the metric of the noncircular (the most general) stationary axisymmetric spacetime in a covariant 
fashion, we adopt the (2 + 1) + 1 formalism developed by Gourgoulhon & Bonazzola (1993). Note that this formalism 
is different from the (2 + 1) + 1 formalism by Maeda, Sasaki, Nakamura, & Miyama (1980) and by Sasaki (1984; see 
also Nakamura, Oohara, & Kojima 1987), which is suitable to the axisymmetric gravitational collapse. Here we adopt 
the formalism by Gourgoulhon & Bonazzola (1993) because it is more convenient for a spacetime which is not only 
axisymmetric but also stationary. 

Let n M be the unit timelike 4- vector orthogonal to the t — const, hypersurface T, t and oriented in the direction of 
increasing t, 

n„ = -Nt tlx . (21) 

The 3-metric induced by g^ v on S t is given by 

h^ v = g^ v + n M n„. (22) 

Similarly, let be the unit spacelike 4- vector orthogonal to the t — const, and ip — const, hypersurface T, tip and 
oriented in the direction of increasing ip, 

= Mhjtp^ = Mip^, (23) 
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where the vertical stroke | denotes the covariant derivative associated with the 3-metric /i M „. The induced 2-metric on 
St p is given by 

H^u = - m^rriv = g^ u + n^n v - m tl m v . (24) 

The covariant derivative associated with the 2-metric is denoted by a double vertical stroke ||. There is a relation 
between the determinants as ^/^g = N^h = NM^/H. 

Any 4- vector can be decomposed into its projection onto T, tif> , the component parallel to and that to m^. The 
Killing vectors are decomposed as 

T]^ = Nn^ -N» = Nn^ - MN v m^ - , (25) 
f = Mm' 1 - Me", (26) 

where the shift vector is (minus) the projection of 77^ onto St, Ms M is (minus) the projection of ^ onto £ t¥ ,, and 
is the projection of onto T, ttp . For our choice of the coordinates, i.e., for x° = t and x 3 — (p, the component 
expressions for n M and are 

1 N 1 N 2 N v \ 

N'~N'~N , ~N~J , (27) 
Ms 1 A/ s 2 1 \ 

°'1w"'1w"'mJ- (28) 

Note that = (0, iVs 1 , iV s 2 , 0) and = A?£ A + N v Ms A . We can express the 4-metric g^ v in terms of N, N v , 
Nx A , M, M S A and H AB as 

g^dx^dx" = - [iV 2 - M 2 (A^) 2 - N SA N^ A * 2 - 2 (Af 2 A^ - N^ A M^ A ) dtdip 

- 2N BA dtdx A + H AB dx A dx B - 2Mv A dipdx A + (A/ 2 + M E aM e a ) c^ 2 , (29) 

where the functions 2V, A^, A^" 4 , M, M^, A and i^s depend only on the coordinate (x 1 ,^ 2 ). We note that the 
presence of Ae m and characterizes the degree of noncircularity of the spacetime. Since we only assume that 

physical quantities are independent of x° = t and x 3 — ip, the metric <? M „ in equation (|29|) has some degrees of freedom 
in the choice of coordinates. Nevertheless, the norms N^N^^ and Me^Ms^ cannot be set equal to zero for a general 
noncircular spacetime. In § 12.51 the GS equation will be given as an equation projected onto T, tv . 

2.4. Physical quantities from flux function ^ 

Provided that the metric g^ is given and the conserved quantities E(^) (or D(^)), L(*ff), C(\&) and S(^f) are 

given as functions of 'J, all the physical quantities can be evaluated once the (effectively 2-dimensional) configuration 
of the flux function \& is known (Ioka & Sasaki 2003). 

The fluid 4- velocity is expressed as 

u» = u,(t^ + fie) + u 6 (e + Gr)») + u^, (30) 



where 



E-Q.L _ D 

Gr/fJ- Grj Grjl-l G r j 

[ ' M s 

G 6l i \c n C 2 P ) [ l G n C 2 P ) + G 



(31) 



and we define 



(L - QE) ( 47T/X \ / 47T/X V . . 
M C = f^— (7^27)1 !- 7^-^27) +— ' 

e Uv' + W) = .903 + »ff33 m) 

r]^ + Q^) .9oo + ttg 03 ' 

G v = -{^ + + ^) = -(.9oo + 2Qg 03 + tt 2 S 33 ), (35) 

- (£„ + e VlM )(e + ®if) = 333 + 29 ff o3 + e 2 ff0 o, (36) 

Nx = u^(N^ + nM^), (37) 

Mv=u^{M^ + QN^), (38) 

= e^ afs n a m . (39) 

In equation (fTTTTf) - the orthogonality {rf + fi^X^ + = is satisfied. Note that Mau 2 ■= 4Trp,/G v C 2 p is the 

square of the effective Alfven Mach number Mau- At the Alfven point Ma\i = 1, the numerator L — QE in equation 
1321) should vanish to keep the velocity finite. 
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Thus, from equations l|30(l - (|33|l . given the metric g^ and the conserved functions E(*f?) (or D(^f)), L(^f), 0($), 
C(\&) and <S(^), if the density p and the enthalpy p are additionally known (see below), the fluid 4-velocity can 
be obtained from the flux function iff and its first derivatives *ff,A- The magnetic field is also calculated from the flux 
function Hf by using equation (|15|l . The (2 + 1) + 1 decomposition of the fluid 4-velocity is performed with equations 
(25), (23 and I® as 



where 



u n = N(u v + Qu^), 

u m = M [(ft - N v ) u v + (1 - N v @) ut] , 



(40) 

(41) 
(42) 
(43) 



The pressure p, the internal energy e, the enthalpy p and the temperature T are functions of the density p and the 
entropy S from equations (|19p. (SJ and (5). Hence, given S 1 as a function of 'J', the only remaining quantity to be 
known is the density p. The density p is determined by the normalization of the 4-velocity, which yields 



D 2 



G v p 2 



{Air) 2 {L-QE) 2 
G^CV 



1 



47T/i 

~G^C^~P~ 



N 2 M 2 C 2 p 2 + ~G V 



This equation is called the wind equation. 
The following components of the electric current are also calculated from the flux function, 



1 



J = ■ 



AirNMv H 
1 



NMVHF 



1 



AnNMVH 



NMVHF 



? 3A 



AttNM 
1 

AnNM 



NMF 



0A\ 



>\\A> 



NMF 



3A\ 



'\\A ■ 



(44) 

(45) 
(46) 



where 



g 0B g A0 ) B 



9 0B 9 A3 



^03 „AB 



N 2 

M A N^ B M§e BC 



N 2 M 




jv^B + ^g g -g g ) r 12 



n v h ab _ M A N S B 



M 2 



F\2 

7S' 

2 



(47) 



31„A2 



9 32 9 A1 )F 12 



/ 1 /jyv^ 
I M 2 _ [~N y 



ab Ny A Ny B 
N 2 M 2 



N^N^ B M§ e BC 
N J l " - X' 2 f ' N 2 M 2 

F 12 =C v ^gp(u 3 - flu ) = Cy/^gpu/: (1 - flS) . 

Thus, J° and J 3 are expressed in terms of 'ff and its first and second derivatives. 

2.5. Relativistic Grad-Shafranov equation 
The equation for the flux function "J, that is, the GS equation is given by (Ioka & Sasaki 2003) 

1 



F12 

Vh' 



(48) 
(49) 



j 3 - n j° 



AB 



NMC 



(^usa)\\b - p(u n + Sue) 



dE d(CSl) 



p(v,£ + Qu n ) 



dh A dc 



,rf = o,(50) 



where the double vertical stroke |j denotes covariant differentiation with respect to the 2-metric Hab, and we have 
defined the auxiliary quantity 



A = — (u B 3 - u 3 B ) = -^Cp ( G(.ui: - M s ) (.g 00 + ^303) 



(51) 



In the previous subsection, we have seen that u v , u^, us A , Q, p, p, T, J° and J 3 are all expressed in terms of iff and its 
derivatives, given the conserved functions E(iff) (or D(iff)), L(iff), , C( V E') and S(*ff), and the metric g^ v . Thus, 
in order to obtain the configuration of matter and electromagnetic fields in a curved spacetime, we have only to solve 
the GS equation (|50|l with the aid of the wind equation l|44|). 
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3. GRAD-SHAFRANOV EQUATION IN THE WEAK MAGNETIC FIELD LIMIT 

It is formidable to solve the GS equation (|50|) in a general case because of its high non-linearity. However, it becomes 
tractable if we take the weak magnetic field limit. In this section, we derive the GS equation in the weak field limit 
with the flux function \& being the perturbation parameter. This perturbation method is essentially similar to that 
developed by Chandrasekhar (1933; Chandrasekhar & Miller 1974) and Hartle (1967; Hartle & Thorne 1968) for slowly 
rotating stars, in which the perturbation parameter is the angular velocity. 

3.1. Zeroth order 

In Ioka & Sasaki (2003), we showed that no magnetic field limit is obtained by letting \P — > and C — > oo. In this 
limit, the GS equation H50|) and the wind equation (|44|l are reduced to the integrability condition and the Bcrnouilli's 
equation for the rotating isentropic (dS/d^ = 0) star, respectively (Ioka & Sasaki 2003). Thus a rotating fluid star is 
the zeroth order configuration in the weak magnetic field approximation. 

Let us further suppose that there is no rotation £1 = at the zeroth order in Then the background metric g^y is 
spherically symmetric, 

g^dx^dx" = -e 2v dt 2 + e 2X dr 2 + r 2 (d9 2 + sin 2 9dtp 2 ) , (52) 
and the zeroth order configuration is obtained by the Tolman-Oppenheimer-Volkoff equations, 

™' = 47rr 2 p(l + e), (53) 

p' = - \e 2X pp (to + 47rpr 3 ) , (54) 

v' = — -p', (55) 
pfj, 

where the prime denotes differentiation with respect to r, and the mass m(r) inside radius r is defined by 

e -=(l-^)~\ (56) 

The value r = at which p = is the radius of the star, and the total mass of the star is given by m(R*) = M*. 
The boundary conditions are m(r = 0) = and 2u(r = i?*) = ln(l — 2M*/R*). 

3.2. First order 

We perform linearization with respect to the flux function 4'. To proceed, we have to specify the functional forms 
of the conserved functions D(^f) (or E(^f)), L($>), f2(W), C(4 r ) and S(^), which characterize the configuration of the 
fluid flow and the electromagnetic field. In this paper, we assume the following. 

• C($): We consider the case in which the magnetic field is confined in the interior of a star. Hence, if? must 
vanish at the stellar surface. Then, the Alfven point M| lf = Air p/G v C 2 p = 1 exists near the stellar surface p ~ 
in equation (|32l) unless C diverges at the stellar surface. Therefore, we assume the form 

C=%, (57) 

where C — const. = 0(1). We note that as long as the magnitude of C is much greater than ^?, the discussion 
below is valid, but we assume it to be O(l) for the sake of order counting. 

• fi(^): Since our primary interest is the effect of magnetic fields on the stellar structure, we consider a slow (rigid) 
rotation which is of the same order as the one induced by the presence of a magnetic field. Hence we set 

Q = const. = 0(* 2 ). (58) 

Here, although we may consider the case f2 = O(W), we choose not to do so, and assume that the rotation is 
subdominant. Note that the observed magnetars have long rotation periods (~ sec), so that this assumption is 
justified. 

• S(ty): We consider an isentropic star for simplicity, 

S = const.. (59) 

• D(*if) and L(^): These functions are chosen in such a way that the GS equation becomes separable with respect 
to r and 8. Namely, 

D = D a + D^ + L> 2 , (60) 



L = const. = O (C) , (61) 

where Do — const. = O(l), D\ = const. = 0(^) and D2 = const. = 0(^> 2 ). Note that the conserved function 
L(^>) does not necessarily vanish at zeroth order at which there is no rotation. This is because L(^) does not 
exactly coincide with the angular momentum. 
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With these assumptions for the conserved quantities, we can estimate the orders of the variables from equations l|31[l 
Pty. and {G5 as 

9 = 0(* 2 ), u n = 0(l), us = 0(* 2 ), u^~u^ = 0(^ 2 ), B Q = 0(* 3 ), Bi = 0{*), (62) 

where we have used the fact that the metric perturbations are O ( V T /2 ), i.e., in equation 129fl we have 

N v = O (* 2 ) , = O (* 2 ) , M E " = O (* 2 ) , (63) 

and <7o3 = O (^ 2 )j which will be confirmed in §0] Note that N = e v ', M = rsinf?, H rr = e 2A and Hgg — r 2 at zeroth 
order from equations (|29|l and Q52[l. 
Now let us derive the GS equation in the weak magnetic field limit. To 0(^>), equation <|5(J|I is reduced to 

, dE k dC 

J pU ^- pU Z A M=°> (64) 
where, from equations (|46|l and ^48(1. the electric current is given by 

Up to O (* 2 ), equations l(3"5)l. l|3"6"|) . l|3"T|) and ({3"2*)l are reduced, respectively, to 

G v = ~9oo, (66) 
G C =533, (67) 
D 

(-S'oo)m' 



Hence the wind equation (|44|) is reduced to 



47rL , , 

9009^ P 



°~ ° , = (-SboK = l, (70) 



e A ^ 



(-ffoo)^ 2 M 

where we note that goo and 333 include perturbations up to O^ 2 ). Then, at O the second term in equation 
is given by 

dE d(D + QL) _ v . , 

~ PMr 'd* = = ~ pu ^ 1 = 1)1 ' ( 71 ) 

where we have used u v = 1/V — 5oo = e"" at zeroth order that follows from equation Q70JI. Note that we have Do = p.e v 
at zeroth order. Similarly, from equations l|51|l. I|57|) and (|32l) . the third term in equation (|64J) is approximated as 

- pwA— = =- — ^- e- 2v ^>. (72) 

P 4 d«f Anr 2 sin 2 6\Cj 

Then, by combining equations iffiljl. JHU, f7TJ> and ijTSjl. the GS equation valid to O('I') is given by 

p v+\ / 1 \1 /AttT\ 2 

d r (e v - x d r V) + smOdg -dgV + ^ e 2 ^-"^ + 4irr 2 sin 2 6pe 2X - v D l = 0, (73) 

v ' r l \smB ) \ \ C ) 

where we have used the form of the background metric g^ v given in equation l|52|l . Note that, under the assumptions 
for the conserved quantities we have made, the GS equation at this order is an inhomogeneous linear differential 
equation for '5 with the last term proportional to D\ as a source. 

Let us separate the angular variables in the GS equation. We expand the flux function \t by the vector spherical 
harmonics (Regge & Wheeler 1957; Zerilli 1970) as 

^ = Y,Mr)sm9 n (74) 
1=1 

where Pi is the Legendre polynomial of degree i, and we have discarded the ^-dependence of the harmonics because 
of the axisymmetry. Substituting this form into equation l|73|l , we find the equation for the dipole ipi as 

ft + {v 1 - A') ft x + (l 2 e~ 2v - e 2A Vi - 47rr 2 pe 2A -^ 1 = , (75) 

and for the higher multipoles as 



ft>+(v' -\')ft t + (^L 2 e- 2 » ~ £ -^-^je 2X ^ = (£>2), (76) 
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where we have introduced 

Equations 175(1 and 1(76(1 are the basic equations for the matter and electromagnetic fields for the conserved quantities 
of the assumed forms in the weak magnetic field limit. 

Since we are interested in magnetic fields confined in the stellar interior, or interior magnetic fields with magnitude 
much larger than exterior magnetic fields, we assume that the magnetic field vanishes at the stellar surface. Thus, to 
the lowest order in we require 

MR*) = o, MR*) = o {£ > i). (78) 

Near the origin r — 0, the regular solutions of equations 1(75(1 and ((76(1 should have the behavior 

iP e ^r e+1 a e {£>!), (79) 

where at is a constant. 

First consider the dipole ipi. With the boundary conditions given above, equation l(75() becomes an eigenvalue 
equation with L being the eigenvalue to be determined. There will be a discrete set of eigenvalues. For a given 
eigenvalue, the constant D\ then just determines the normalization of ip±. Now we consider the higher multipoles 
ipt (^ > 2). Once L is determined, it is in general impossible to find non-trivial solutions for £ > 2 that also satisfy 
the same boundary conditions ((78(1 and ((79(1 . Thus, under our assumptions that the conserved quantities are of the 
forms 1(57(1 - 1(61(1 and that magnetic fields are confined in the stellar interior, the only non-trivial solutions are of 
dipole type. Of course, once we relax the assumption of confined magnetic fields, the existence of higher multipoles 
will be allowed. 

In the numerical analysis, it is convenient to integrate equation ((75(1 from the stellar surface with equation 1(78(1 as 
initial data. By varying L, we look for a solution that behaves as regular as possible at the origin. This procedure 
determines the eigenvalue L. Next we integrate equation 1(75(1 from the origin with the behavior ((79(1 and the fixed 
eigenvalue L. Then the constant ai is determined by demanding that the boundary condition at the stellar surface, 
equation 1(78(1 . is satisfied. 

4. METRIC PERTURBATION EQUATIONS 
In this section, we derive the field equations for the metric perturbation, 

Ag^ = g^y - . (80) 

By linearizing the Einstein equations about the background metric g^ 7 we have 

{A 9tnJ , a ' a - A 9tia .y a - Ag^.y 01 + Ag a a ,^) 

+g fll/ {Ag a p' a ^ - Ag a a ,^ - Ag a pR a P) + Ag^R = -16^AT M „, (81) 

where AT^ is the perturbation of the energy momentum tensor, and i?^„ and R are the Ricci tensor and the Ricci 
scalar, respectively. As shown below, the metric perturbation Ag^ is O ( V I' 2 ) • Our method is similar to that developed 
by Hartle (1967; Hartle & Thorne 1968; Chandrasekhar & Miller 1974) for slowly rotating relativistic stars. We obtain 
the field equations for the interior of the star in § 14.11 and solve the exterior equations to give the junction conditions 
at the stellar surface in § 14.21 

4.1. Interior field equations 

The study of the metric perturbation in a spherical spacetime was initiated by Regge & Wheeler (1957) and Zerilli 
(1970). Because of the spherical symmetry of the background, the perturbation equations are separable in the angular 
variables by using ten tensor harmonics. We can eliminate some components of the metric perturbation with the aid 
of the gauge freedom. In the gauge used by Regge & Wheeler (1957), the metric perturbation belonging to a given I 
is of the form, 



f-2e 2u H e It \ 

h fe 4A M £ 
r 2r 2 K e 

\ 2r 2 sin 2 6K t ) 



Pi (cos ( 



/ 7/ 
OWe 1 . a d 

] sm9—P e (cos9), (82) 
\V E W e 



where the first term is of parity (— 1) (even parity), the second term is of parity (— (odd parity) and we may 
discard the (^-dependence of the harmonics because of the axisymmetry. Since there are fewer independent tensor 
harmonics for £ = and £ = 1, we may further simplify the metric by adopting a gauge in which Iq = 0, Kq = 0, 
Vo = 0, W = 0, Ki = and Wi = (Zerilli 1970). 
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The energy momentum tensor perturbation AT),,, is obtained by taking the perturbation of equation (J7|). The 
components of the magnetic field to 0(5') are given from equations (|15|l . (|30|l - Ij33|l . H62(l . 1|63|1 and (I7t)|l as 

S° = 0(* 3 ), (83) 

5 1 = -^-e- A ^iCos6', (84) 

5 2 = -^e-ViSin6>, (85) 

B 3 = ^e->i. (86) 

Since the energy momentum tensor is quadratic in B^, we need AT^ to the accuracy of O {^ 2 )- To this order, from 
equations - E2J>, (ESI and (JTOJ), the fluid flow is given by 

u° = u n = -=t=, (87) 
V~.9oo 



u 1 = -J— er x - v (t/;i) 2 sin 2 6 cos 0, 



7— 



u^-J—e-^^sm'e, (89) 
u 3 = nu v + U£ = Vie-" ~ e~ 2l/ (?/>i) 2 sin 2 9. (90) 



Cp 

Expanding the perturbation of the enthalpy in equation (|19f) in the form 

Ap 



= £ ^ In Kr)hPt (cos 0), (91) 
we obtain from equations l|7U|l and (|82|l the equations 

(Aln M ) = -Ho-|^Vi + ; ^, (92) 

(Aln^) 2 = -ff 2 + ^Vi, (93) 

6 Do 

and (Aln/i)^ = —Hi for £ 7^ 0,2. From the first law of thermodynamics, equation JHJ, with equation (|19|l . we have 
the relations (Ap)e = p(Ap)t and (Ae)^ = (p/ p 2 )(Ap)i, where (Ap)g, (Ae)e and {Ap)i are the £-th pole perturbations 
of the pressure, the internal energy and the rest mass density, respectively. For the polytropic equation of state 
p = Kp 1+1 / n , we have (Ap) e = (1 + l/n)(p/p)(Ap) e . 

From the Einstein equations (|81|l with Ag^ given by equation i|82|) , we obtain the following independent equations 
which have nonzero source terms, 

2 2 n dp . . . 1 / J. 2 -2v 1 2 A /„,, \2 . l-2A/,,/\2 



= 47rrV^(Aln M )o + - + + -e~ 2A (^i) , (Goo) (94) 
H - = -[(Alnp) ]'--?-e-»D 1 ip' 1 

= ^ + 8^e 4A Af + 47rre 2 V(Aln A1 )o + ^ e 2A (^P e - 2l '-^(^i) 2 + ^(V' 1 ) 2 , (G u ) (95) 

Ji = _3|r| eA(V)i)2) (Goi) (96) 

/ 3 = ^t e A (^) 2 , (Gbi) (97) 

Vi" - (1/ + A') V{ - - (V + X' + -) V, = ^Lipe 2X -» (Vi) 2 - 16 W r 2 e 2A tt, (G 03 ) (98) 

V? - W + A') Vi - I (V + X' + \ + ^e 2A ) V 3 = -^-ipe^ {^f , (G 03 ) (99) 

Wa = -|ie A -"(Vi) a , (G 13 ) (100) 

H 2 + i^Ma = |e- 2X - |i a e- 2v (W 2 , (G 22 - G 33 / sin 2 fl) (101) 
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K' 2 - ^e 2X K 2 + -H' 2 - ^e 2X H 2 - - , , 
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1 



87175 



e 4A M 2 



1 



47re^ A pM(Aln/i) 2 - r^e 
3r z 



2A 























Vr 





„2A , 



(G11) 



3r 2 



(102) 
(103) 



where we have indicated the component of the held equation in the parentheses at the end of each equation from 
which it is derived. The remaining components of the Einstein equations that are not listed in the above are either 
trivially satisfied or not independent from the above set because of the Bianchi identities (Zerilli 1970). To avoid the 
cancellation of significant digits in the numerical calculation, it is convenient to define 



Y 2 = H 2 + K 2 V 2A W) 2 - |-e- 2 - ^ (V>i) 2 . 

or 3r* 



(104) 



and solve for Y 2 instead of K 2 . With equations l|93|) and i|101|) . we can rewrite equations l|102[) and (|103(l as 



Y{ + 2v'H 2 



-2A 



«4) 



-2v 



+ -(1/' + A') e- 
r 



-2A 



4to 



H' 2 



2Ai 



2v' 



1 



2 A 



-2A 



MY 



rV'iVh 



4:TTpil 
1 



2m 



Ho 



3 



-2v' 



,2A 



r2 -2l/ 

L e 



/ / \2 8tt 



2A-y 



£>i^i. 



(105) 



(106) 



We have six differential equations ||SSJ|, ffigjl. (|9"§|) . (|105|l and I|1U6|) for six unknown functions; M , (Aln/Lt) , Vi, 
V3, Y 2 and if 2 . A practical method to solve these equations is discussed in the next section. The functions Ho, 1%, I3, 
W 2 , M 2 , K 2 and (Aln/x) 2 are determined algebraically from equations ®, 113, C^UI) - (fTCTTf t . (|TUl|l and (|9"5j> . 

respectively. 

4.2. Exterior solutions and boundary conditions 

The equations derived in the previous subsection apply to the interior of the star. Outside the star, in the vacuum, 
we have ipi = 0, m(r) — M*, e 2l/ = 1 — 2M*/r and Aln/i = and the resulting equations can be solved explicitly. 
The solutions for which the metric g^ v is asymptotically flat are obtained as 

: AM*, 
AM, 



M 
H 



r - 2M* 

h=h = o, 

2AJ 



Vi: 

V 3 - 

W 2 - 

H 2 



r 
AV 

: 8M| 

= 0, 



(— 

\2M* 



5AQ 

8M3 ^ 2 



^(mT- 1 



V 2 = - 



M 2 = - 



5AQ 



2M, 



8M| [r(r-2M»)] 1 / 2 
(r-2M t )ff 2 , 



r 

mT 



where AM*, A J, AV and AQ are constants, V/ is the function defined by 



Vf(z) 



5 + 30z - 210z 2 + 180z 3 + 60z 2 (3z 2 - 5z + 2) In ( 1 - 



and Q™ are the associated Legendre functions of the second kind 

z (5 



Qa(*) = 

QI(Z): 



- 

2-3z 2 



3z 2 ) 3 z + 1 

+ - [z 1 - 1) In 

■ 1 2 1 



(z 2 



a/2 



2 
,1/2 



In 



1' 
z + 1 
z- 1' 



(107) 
(108) 
(109) 
(110) 

(111) 
(112) 
(113) 

(114) 
(115) 

(116) 

(117) 
(118) 
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In the limit z — > oo, we have Vf(z) -> 1/z 3 , Q%(z) -> 8/5z 3 and Q\{z) -> -2/5z 3 . We may identify AM* and AJ 
with the mass shift and the angular momentum of the star, respectively, while AQ and AV^ with the mass quadrupole 
and current hexapole moments, respectively (Thorne 1980), apart from numerical factors. 

Let us now explain our strategy to solve six unknown functions Mo, (Aln/i)o, Vi, V3, Y 2 and H%. First, to obtain 
Mo and (Aln/x)o, we need to specify what quantity is to be fixed when we add the perturbation. Although we may 
fix the central density (Hartle 1967; Hartle & Thorne 1968; Chandrasekhar & Miller 1974) or the total mass, here we 
adopt the baryon mass as the fixed parameter. From equation the total baryon mass is given by (Hartle 1967) 



-gpu 



Then, from equations l|82[) and JS7J, the perturbation of the total baryon mass up to 0{^ 2 ) is given by 

AMb = I ar^irr-e' 



-e 2X P M + (Ap) 
r 



This may be written in the differential form, 

AM' B (r) = 47rrV 



,.,2 A 



pM 



dp 

pp.— Aln/i)o 
dp 



(119) 



(120) 



(121) 



where AMs(r) is the baryon mass within radius r, and we have used the equality (Ap)o = pp(dp/dp)(A ln/x)o for 
an isentropic fluid. To fix the baryon mass, we solve this equation with the boundary conditions AMg(O) = and 
AMb{R*) = 0, together with equations l|94|l for Mo and (|95|l for (Aln/i)o. Near the origin these functions have the 
behaviors, 

1 J, 9 dfc 
dp c 



Mo 



Amp^pc 2 ^- (Aln/i)oc + 2a 2 



(Aln/i) — >(Aln p) 0c , 

AM B ^^-r 3 p cPc ^{A\np) 0c , 
3 dp c 



(122) 
(123) 
(124) 



where ol\ is the constant introduced in equation l|79l) . i.e., ipi ~ * r 2 ot\, and the subscript 'c' means the value at the 
center of the star. With these initial behaviors, equations l|94|l . 195|l and (|121fl can be integrated from the center to the 
stellar surface. We first determine the central value (Aln/i)o c so as to satisfy the boundary condition AMb{R*) = 0. 
Then the value of Mo at r = i?* gives the mass shift AM* from equation l|107|) . and the solution for (A In p)o, together 
with equations (|92|l and (|108fl . determines the constant D^. 
Next, to obtain V%, it is convenient to express it as the sum of a particular solution and a homogeneous solution 



- C1 V} H \ (125) 

where c\ is a constant and the homogeneous solution V^ H ^ satisfies equation l|98|l with the right-hand side being zero. 
Near the origin, we have 

V { P) ^r 2 , v[ H) ^r 2 . (126) 

Then, V} P) and v} H) can be obtained by integrating equation (|98|l from the center with the initial behavior (|126fl 
to the stellar surface. The constant c\ in equation l|125|) and the angular momentum AJ in equation (|110f> can be 
determined by requiring the continuities of the value V\ and its derivative V{ with the exterior solution in equation 

rrm . 

A similar method can be applied to equation 1)99(1 for V3. In this case, the behavior near the origin is 

V 3 {P) ^r\ v( H) ->r\ (127) 

Then the solution determines the current hexapole moment AV from equation (|lll|l . 
Finally, to obtain Y% and Hi , we again express them as the sum of a particular solution and a homogeneous solution 



Y 2 = Y 2 



c 2 yI h \ H 2 =Hf>+c 2 H^, (128) 

where c 2 is a constant and the homogeneous solutions and satisfy equations 1|105|) and (|106fl . respectively, 

with the right-hand side being zero. Near the origin, they behave as 

„4 



(P) 



AH) 



Y 



(P) 



Y 



(H) 



6 V 

• — t t \p 



2 — 1v r 2 o 

°a 1 — 8np c e 



-Dion 



H. 



(P) 



"22 

3 u 



p c e c + 3p c ) , H. 



(H) 



(129) 
(130) 



These functions can be calculated by integrating equations l|105fl and 1)1 06f) from the center with the initial behav- 
iors (|129|) and (|130fl to the stellar surface. The constant c 2 in equation l|128fl and the mass quadrupole moment AQ in 
equa tions l|113|) and l|114|) can be determined by requiring the continuities of Y 2 and H 2 with the exterior solutions (|113H 
and ifHifl . 
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5. NUMERICAL RESULTS 

In this section we present the results of our numerical integration of the linearized GS equation (|75|l and the metric 
perturbation equations in § 0] We adopt the fifth-order Runge-Kutta method with the adaptive stepsize control (Press 
et al. 1992) and use it iteratively when necessary. We adopt the polytropic equation of state 

p = Kp 1+1 ' n , (131) 
with the polytropic index n = 1. We consider the non-rotating case f2 = to concentrate on the magnetic effects. 

5.1. Magnetic field and fluid flow structure 
The GS equation (|75|l for ipi is the basic equation to obtain the fluid flow and the electromagnetic field. The 
boundary condition 178fl is satisfied for discrete eigenvalues L. The first few roots of the (dimensionless) eigenvalue 
are given in TableQfor several relativistic factors M*/i?*. Although L can take both signs, we assume L > in 
the following. 

The components of the magnetic field B M are given by equations - ijSfill . Although the magnetic field is defined 
in the fluid frame in equation (0, the fluid 4- velocity there may be identified with the unit vector n M normal to 
the t =const. hypersurface to the accuracy of 0(\&). Therefore, it is convenient to express the magnetic field in terms 
of the orthonormal tetrad components in the background metric (I52f) . They are given by 

2 

5 f = -— V>iCOs6>, (132) 
B^-e'^smO, (133) 

B = -e- v ^) 1 sm.e. (134) 
r 

In Figures n we show Bf(9 — 0), Bg(9 = tt/2) and B^(9 = ir/2) normalized by 2a\ as a function of radius r for the 
smallest two eigenvalues of L. Note that the magnetic field strength at the center of the star is given by B 2 — (2ai) 2 . 
From these figures, we can see that the magnetic structure is more complicated for a larger eigenvalue L. To be precise, 
the number of nodes at which the r-derivative of the magnetic field vanishes increases as L becomes large. This is also 
illustrated in Figures [21 in which the magnetic field lines projected on the meridional plane, i.e., the W =const. lines 
in the (r, #)-plane, are shown. In Figure [21 we show the projection of a truncated piece of a magnetic field line onto 
the equatorial plane. Only the part of it in the upper hemisphere (9 < tt/2) is shown. As we can see from the figure, 
the flux surface forms a torus and the magnetic field line winds around the torus. 

The fluid flow is determined from equations (|87J) - l|9Q(l . The orthonormal tetrad components of the fluid velocities 
in the background metric l|52|) are given by 

2 ""{ipif sin 2 9 cos 9, (135) 



Cpr 2 
1 



vx = - — 



Cp 



r 



■e- x - v M' 1 sin 3 9, (136) 



v = n e -' , rsm9- -^e~ 2v (V>i) sin 3 9, (137) 
Cpr 

where we set f2 = in the present analysis. In Figures 0] we show v?{9 — 7r/4), vg(9 = tt/2) and v^(9 = tt/2) as 
functions of radius r for the smallest two eigenvalues of L, normalized by the dimensionless quantity, 

M^^|_ (138) 
C Ml 

As in the case of the magnetic field, the velocity field is more complicated for a larger eigenvalue. In Figure[5l we plot 
the velocity field on the meridional plane (ip =const. plane). The flow line resides on the constant flux surface. 
In equation i|138|l . the dimensionless quantity, 

tj . RWx 3 (g c 2 /87r)(47ri?g/3) 

describes the ratio of the magnetic energy to the gravitational energy, apart from a numerical factor. The quantity C 
controls the magnitude of the meridional flow. As mentioned in § 13.21 our perturbation method is valid as long as 

|C|»|i? 3 a 1 | = 0(vI/). (140) 
The limit \C\ — > oo corresponds to the limit of zero meridional flow. 
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5.2. Mass shift and magnetic helicity 

By integrating the metric perturbation equations in § ^ according to the prescription of § 14.21 we can calculate the 
change in the central density (Ap)o c , the mass shift AM*, the angular momentum A J, the mass quadrupole moment 
AQ and the current hexapole moment AV^. In Tabled we list the numerical results for several configurations. Each 
value is normalized appropriately as indicated in the first column. Here the dimensionless variable TZy defined in 
equation (|138f) measures the magnitude of the meridional flow velocity and TZm defined in equation 1|139|) measures 
the ratio of the magnetic energy to the gravitational energy. We also introduce the dimensionless magnetic helicity 
H.M as defined in equation 1)145(1 below. The numbers given in Table can be easily converted to physical units once 
we specify the mass M* and the dimensionless variables TZy (or C) and TZm (or a{). 

The mass shift AM*/M*Hm increases for larger eigenvalue L. This means that the configuration with the smallest 
eigenvalue L = L\ is the lowest energy state for a fixed baryon mass Mb and total magnetic helicity TLm- The magnetic 
helicity is conserved in a perfect conducting plasma (Woltjer 1958), and it describes topological properties such as 
links, twists and kinks of magnetic field lines (Moffatt 1969; Berger & Field 1984). When the topology of the field lines 
changes, the conductivity should be finite, which implies non-conservation of the magnetic helicity. However, it was 
conjectured that the total magnetic helicity is still approximately conserved on reconnection or relaxation time scales 
(Taylor 1974; Berger 1984). This hypothesis has been successful in laboratory plasmas (Taylor 1986). It is uncertain 
whether or not the total magnetic helicity is nearly conserved in the neutron star. In any case, we may speculate that 
the configuration with the smallest eigenvalue is preferred if the free energy could be released by the rearrangement of 
the magnetic field (see also Ioka 2001). 

In the covariant language, the 4-current of magnetic helicity is defined (Carter 1978; Bekenstein 1987) by 

= \e^A v F a p. (141) 

This is conserved in the ideal MHD, 

= 0. (142) 

Thus the total magnetic helicity is defined by J d 3 Xy/ —gH° . Although the 4-current of magnetic helicity is not 
gauge-invariant locally, the total magnetic helicity is gauge-invariant for confined magnetic fields. From equations (|14ll . 
JH3 and (JUl we have 

F 12 = Le x ~>isin0. (143) 
Hence, using the gauge freedom to take A 2 = 0, we find 

A x = Le A ~^icos6>. (144) 
Therefore, we define the (dimensionless) total magnetic helicity by 

H M ■■-JpJ d 3 x^H a = ± 1^ dr^Le x - (^ f . (145) 

5.3. Spherical and quadrupolar deformations 

The magnetic stress distorts a star. The isobaric (p=const.) surface that lies at a radius ro in the background 
configuration is displaced in the perturbed configuration to a radius 

r = r + (Ar)o(ro) + (Ar) 2 (r )P 2 (cos#), (146) 

where, from equation l|55ll and (Ap)» = p(Afi)^, we have 

(Ar) e = i-(Aln^ (^ = 0,2). (147) 
v 

The spherical deformation is characterized by the monopole displacement (Ar) . However, to describe the quadrupolar 
deformation in a geometrically invariant manner, the coordinate radius r is not quite appropriate. Instead of r, the 
circumferential radius r 2 — r 2 [I + 2K2P2 (cos 9)} is more appropriate (Hartle & Thorne 1968). Then the isobaric surface 
is expressed as 

r*(r , 9)=r + (Ar) (r ) + [(Ar) 2 (r ) + r K 2 (r )} P 2 (cos6>). (148) 
We describe the quadrupolar deformation by the ellipticity defined by 

"(Aln^) 2 



r„(r,7r/2)-r„(r,0) 3 
<(r) = 



(149) 



r 2 

In Table [3 we also list (Ar)o and e* at the stellar surface r = i?* for various configurations. We can see that the 
stars are more deformed for larger eigenvalues of L, for the baryon mass and total magnetic helicity fixed. This is 
intuitively consistent with the results in the previous subsection that the energy (the mass shift) is higher for larger 

L for constant baryon mass and total magnetic helicity. We also find that the configurations are prolate rather than 
oblate since e* < 0. This is in contrast with the case of stars with only poloidal fields which are oblate rather than 
prolate (Bocquet et al. 1995; Konno, Obata, & Kojima 1999). This difference is due to the fact that the magnetic 
stress of a toroidal field tends to make a star prolate, working like a rubber belt tightening up the equator of the star 
(Ostriker & Gunn 1969; Ioka 2001). Our model is the first example that demonstrates the validity of this picture in 
relativistic stars. 
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5.4. Frame dragging 

In general relativity, inertial frames are dragged not only by the fluid flow but also by the magnetic field. The 
dragging of inertial frames manifests itself in the off-diagonal components of the metric 1|82|) . In addition to the tip- 
component (Vi) which is familiar in the case of rotating stars and Kerr black holes, there also exist the ir-component 
(It) and the rc/?-component (Wt) in the Regge- Wheeler gauge. These components originate from the meridional flow 
and the toroidal magnetic field, which are incorporated into the relativistic stellar model for the first time in this 
paper. 

The 4- velocity of an observer whose world line is orthogonal to the t =const. hypersurface St is given by n M in 
equation J2U- To the lowest order in the coordinate angular velocity and the coordinate radial velocity of such 
observers are given by 

— = -jr = iV 1 = -e~ 2 x h cos 6 - e~ 2A / 3 - (5 cos 3 9 - 3 cos 9) , (151) 
at n" 2 

from equations (2JJ, and (|g2|l . 

In Figures IH1 we show the ^-component of the metric perturbation that describes the coordinate angular velocity, 
for M*/i?* = 0.2 (relativistic case) and M*/i?* = 0.01 (Newtonian case). The left panel shows the dipole (i = 1) 
component of the (normalized) coordinate angular velocity (Vi/r 2 )/(M*TZv / Rl) and the right panel shows the hexapole 
(£ = 3) component (V3/V 2 )/ '(M*1Zv / R*)- From these figures, we can see that the normalized coordinate angular 
velocity does not depend on the relativistic factor M*/i?* so much. Thus, if the velocity of the meridional flow is 
fixed, i.e., IZy =const., the coordinate angular velocity is nearly proportional to (M*/i?*) 2 for fixed mass. Note that 
the i(/?-component (Vt) is gauge-invariant for stationary perturbations (see equation (D3) in Zerilli 1970). 

In Figures we show the dipole and hexapole components of the (normalized) coordinate radial velocity in equation 
11511) . — e~ 2X h / (M*lZv / R*) and —e~ 2X Iz/{M lt 'R,v/R*), respectively. From these figures, and from the left panel of 
Figures 0] where we see that the fluid velocity v is about v ~ O.llZy, as an order of magnitude estimate, the coordinate 
radial velocity is found to be ~ (M*/R*)v. It is interesting to note that the presence of the components I\ and I3 
violates the reflection symmetry about the equatorial plane, i.e., 1% and I3 are of parity (—1) in equation l|82|). Note 
that the ir-component (Ig) is gauge-specific by itself (see equation (D4) in Zerilli 1970). For example, we may choose 
a gauge in which goi — 0. However, for such a gauge we will have 502 7^ 0. What is significant is that we have a new 
frame dragging effect due to non-vanishing iVs M because of the noncircularity of the spacetime, as pointed out at the 
endof§HC. 

In analogy with the above arguments, let us consider the components of the unit spacelike 4-vector m* 1 orthogonal 
to the t =const. and ip =const. hypersurface in equation (|23J) . To the lowest order in 'J, from equations (|28[) . 
(f29l and (JH2J), we have 

^- = ^4 - Mv 1 = 3e~ 2X W 2S m 2 9cos9. (152) 

dip mr 

In Figure|Sl we show the quadrupole component of the (normalized) coordinate derivative e~ 2X Wij '(M 2 TCm / f R*) for 
M*/ii* = 0.2 and M*/i?» = 0.01. We find from this figure that dr/dip oc M^/i?* for fixed mass and total magnetic 
helicity. This 'frame dragging' effect originates from the magnetic field rather than the fluid flow, since the right-hand 
side of equation l|100|) does not depend on C (or TZy)- The component W2 also violates the reflection symmetry 
about the equatorial plane, i.e., W2 is of parity (—1) in equation l|82l) . Note, again, that the r^-component (Wi) is 
gauge-specific (see equation (D3) in Zerilli 1970). For example, we can take a gauge in which 513 = 0. However, once 
again, an important point is that we have non- vanishing because of the noncircularity, which is a gauge-invariant 
statement. 

6. SUMMARY 

We have investigated stationary axisymmetric configurations of magnetized stars in the framework of general rel- 
ativistic ideal MHD. A toroidal magnetic field and meridional flow in addition to a poloidal magnetic field are in- 
corporated into the relativistic stellar model for the first time. We have treated the magnetic field and meridional 
flow as perturbations, but no other approximations were made. We have restricted ourselves to a star with a dipole 
magnetic field, slow rotation and polytropic equation of state, and considered magnetic field configurations confined 
in the interior of the star. 

The solutions are found to be characterized by discrete eigenvalues that describe the radial profiles of the magnetic 
field and fluid flow. Namely, the magnetic field and velocity field configurations become more complicated for larger 
eigenvalues. We have found that the stellar shape is prolate rather than oblate because of the stress of the toroidal 
magnetic field. For fixed total baryonic mass and magnetic helicity, more spherical stars are found to have lower energy. 
We have also found two new types of frame dragging that differ from the one in rotating stars and Kerr black holes. 
These effects are due to the presence of the meridional flow and toroidal magnetic field. Interestingly, these effects 
violate the reflection symmetry about the equatorial plane. This means that the magnetic field or the meridional flow 
naturally selects the preferred direction about the axis of symmetry. This mechanism may be related to the neutron 
star kick or the jet formation. 
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Table 1. Th e fi rst few roots of the (dimensionless) eigenvalue R*L 

IN EQUATION 1751 ARE SHOWN FOR SEVERAL RELATIVISTIC FACTORS M*/R*. 

Note that the configuration M*/R* = 0.24 is unstable against radial 

PERTURBATIONS. 



M*/R, = 0.24 M,/R, = 0.2 M*/R* = 0.1 M*/R* = 0.01 



|i?„,Li 2.8704 3.9079 5.7924 7.2638 

\R*L 2 4.1284 5.6017 8.3151 10.475 

jiJ«£ 3 5.3393 7.2568 10.795 13.615 

\R*L 4 6.5432 8.9031 13.258 16.729 

\R*L 5 7.7448 10.545 15.711 19.829 



Table 2. Parameters of weakly magnetized stars in general relativ- 
ity ARE SHOWN FOR SEVERAL RELATIVISTIC FACTORS M„/i?„ AND EIGENVALUES 

L given in Table Shown are (Ap)o c (central density shift), AAi* 
(mass shift), A J (angular momentum), AQ (mass quadrupole moment), 
A V (current hexapole moment) , (Ar)o (1 = deformation of the stel- 
lar surface), e* (ellipticity at the stellar surface), and Hm (dimen- 
sionless MAGNETIC HE LICITY ) . TH E DIM ENSIONLESS QUANTITIES IZy AND TZm 
DEFINED IN EQUATIONS 11381 AND 11391 . RESPECTIVELY, MEASURE THE MAGNI- 
TUDE OF FLUID FLOW AND MAGNETIC FIELD ENERGY. THE INTEGER SUFFIXED 
TO EACH ENTRY DENOTES THE POWERS OF 10 BY WHICH THE ENTRY NUMBER 
MUST BE MULTIPLIED (E.G., 1.1-2 MEANS 1.1 X 10~ 2 ). 
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Fig. 1. — The components Bf. (8 = 0), Bg(0 = tt/2) and B^(9 = 7r/2) of the magnetic field in the orthonormal basis in equations 11321 
— 11341 normalized by the central magnetic field 2a\ are shown as functions of radius r/R t . We adopt the polytropic index n = 1, the 
rclativistic factor M*/R* = 0.2. The left panel shows the case of the eigenvalue L = L\ and the right panel the case of L = L2 given in 
Table □ 




Fig. 2. — The magnetic field lines projected on the meridional plane, i.e., the \P =const. lines in the (r, 0)-plane. The dashed line is the 
stellar surface. We adopt the polytropic index n = 1, the relativistic factor M*/i?« = 0.2. The left panel shows the case of the eigenvalue 
L = L\ and the right panel the case of L = L2 given in Table HI 
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-1 -0.5 0.5 1 

r sin#cos<^/R* 

Fig. 3. — A truncated piece of the magnetic field line on a certain flux surface =const. surface) projected onto the equatorial plane 
(9 = n/2 plane). Only the part of it in the upper hemisphere (9 < 7r/2) is shown. The dashed line is the stellar surface. The dotted lines are 
the cross section of the flux surface with the equatorial plane. We adopt the polytropic index n = 1, the relativistic factor M*/i?* = 0.2, 
and the eigenvalue L = L\ in Table HI 
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Fig. 4. — The components Vf (9 = 7r/4), Vq(8 = tt/2) and v^{6 = tt/2) of the fluid velocity in the orthonormal basis in equations 11351 — 
11371 . normalized by the dimensionless variable TZy in equation 11381 . as functions of radius r/R*. We adopt the polytropic index n = 1, 
the relativistic factor M-*/R t = 0.2. The left panel shows the case of the eigenvalue L = L\ and the right panel the case of L = L2 given 
in Table 
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Fig. 5. — The velocity field in the m eridio na l plan e [tp =const. plane). The lengths of the vectors are proportional to the fluid velocity 
in the orthonormal basis in equations 11351 - 11371 . The dashed line is the stellar surface. We adopt the polytropic index n = 1, the 
relativistic factor M*/R* = 0.2, and the eigenvalue L = L\ in Table HI 
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Fig. 6. — The tip-component of the dipole (I = 1) metric perturbation (left panel) and that of the hexapole (Z = 3) metric perturbation 
(right panel) as functions of radius, namely, Vi/r 2 normalized by M^R.y I R\ and Vs/r 2 normalized by M*1Zy / R 2 , respectively. The solid 
and dashed lines are for M*/R* = 0.2 (relativistic) and M*/i?« = 0.01 (Newtonian), respectively, while the bold and thin lines for L = L\ 
and L = L2 in Table respectively. We adopt the polytropic index n = 1. These metric components describe the standard frame dragging 
due to the angular momentum. 




Fig. 7. — The tr-component of the dipole (£ = 1) metric perturbation (left panel) and that of the hexapole (I = 3) metric perturbation 
(right panel) as functions of radius, namely, —e~ 2x I\ normalized by M t 1Zy / R* and — e~ 2A i3 normalized by M^IZy respectively. The 
solid and dashed lines are for M*/i?* = 0.2 (relativistic) and M t /R„ = 0.01 (Newtonian), respectively, while the bold and thin lines for 
L = L\ and L = L2 in Table 111 respectively. We adopt the polytropic index n = 1. These metric components describe a new type of frame 
dragging due to meridional flow. 



22 Ioka & Sasaki 




Fig. 8. — The r</j-component of the quadrupole it = 2) metric perturbation as a function of radius, namely, e 2X W2 normalized by 
M^TiM / R*- The solid and dashed lines are for M*/i?» = 0.2 (relativistic) and M*/il* = 0.01 (Newtonian), respectively, while the bold 
and thin lines for L = L\ and L = L2 in Table HI respectively. We adopt the polytropic index n = 1. This metric component describes a 
new type of frame dragging due to a toroidal magnetic field. 



